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Abstract 

In this paper, we study the effects of polynomial f{R) model on 
the stability of homogeneous energy density in self-gravitating spher¬ 
ical stellar object. For this purpose, we construct couple of evolution 
equations which relate the Weyl tensor with matter parameters. We 
explore different factors responsible for density inhomogeneities with 
non-dissipative dust, isotropic as well as anisotropic fluids and dissipa¬ 
tive dust cloud. We find that shear, pressure, dissipative parameters 
and f{R) terms affect the existence of inhomogeneous energy density. 
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1 Introduction 

Recent cosmological evidences predicted by different measurements indicate 
transition of our universe from matter dominated epoch to accelerating ex¬ 
pansion state (Perlmutter et al. 1999, Reiss et al. 2007, Komatsu et al. 
2011). The accelerating cosmic expansion has been prompted by an enig¬ 
matic ingredient with large negative pressure, dubbed as dark energy. To 
explain its nature, different models like cosmological constant, phantom, 
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quintessence, Chaplygin gas etc. have been established. The exploration 
of modihed gravity theories obtained by modifying geometric gravitational 
part of Einstein-Hilbert action has received much attention in mathematical 
physics. The f{R) gravity (Capozziello 2002, Nojiri and Odintsov 2011) is 
one of the most viable theories in which Ricci scalar is replaced by its non¬ 
linear generic function. Among important features of this theory, the likely 
one is to present a model that represents early as well as late-time universe 
expansion in the absence of dark component. Bamba et al. (2012) intro¬ 
duced unified model for inflation as well as late cosmic expansion model in 
this theory. There exist number of f{R) models (Faraoni and Nadeau 2005, 
Nojiri and Odintsov 2007, Hu and Sawicki 2007, Bamba et al. 2012) that 
correspond to cosmological constraints and pass experimental test. 

Anisotropic pressure in matter conhgurations results from several astro- 
physical factors like pion condensation (Sawyer 1972), various types of phase 
transitions (Sokolov 1980), presence of a solid core, superfluids (Heiselberg 
and Jensen 2000) as well as strong magnetic held (Yazadjiev 2012). It is 
worth stressing that for stable huid distribution, anisotropy can endorse out¬ 
wardly increasing density within the core of star (Horvat 2011). After the 
ground work of Bowers and Liang (1974), there have been number of papers 
on pressure anisotropy (Herrera and Santos 1997, Bohmer and Harko 2006, 
Herrera et al. 2008, Sharif and Yousaf 2012, Mimoso et al. 2013, Sharif and 
Bhatti 2013a, 2013b, 2014b) which assert that anisotropy may have non- 
negligible consequences on the structure and properties of self-gravitating 
systems. Thirukkanesh and Ragel (2013) presented spherically symmetric 
compact star models with anisotropic pressure which help to understand 
strange quark stars. 

Gravitational collapse is the phenomenon in which massive body falls 
inward due to the action of its own gravity that may lead to stars, star 
clusters and galaxies from interstellar gas. This occurs due to extremely in¬ 
homogeneous initial state thereby showing the importance of energy density 
inhomogeneities in the collapse process. Penrose (1979) laid much empha¬ 
sis on the importance of energy density inhomogeneity in the gravitational 
time arrow by relating inhomogeneous density with Weyl tensor. Eardley 
and Smarr (1979) asserted that inhomogeneous spherical dust conhguration 
leads to naked singularities for inhomogeneous collapse. Herrera et al. (1998) 
discussed the role of density inhomogeneities and local anisotropy of pres¬ 
sure on the structure and evolution of spherically symmetric adiabatic self- 
gravitating objects through the active gravitational mass. Further, Herrera 
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et al. (2004) investigated density inhomogeneity effects on the evolutionary 
phases of dissipative anisotropic spherical systems by evaluating a link be¬ 
tween the Weyl tensor and local anisotropic pressure. Bamba et al. (2011) 
discussed matter instability and curvature singularity in the star collapse 
with f{R) background. 

Ziaie et al. (2011) studied collapsing mechanism of a star satisfying 
barotropic equation of state in f{R) theories and found hnite-time singu¬ 
larities. Borisov et al. (2012) analyzed spherical collapse in metric f{R) 
gravity with the help of time evolution numerical simulations. Guo et al. 
(2014) studied collapse of spherical star in Einstein f{R) frame and con¬ 
cluded that this may lead to de-Sitter Schwarzschild black hole. We have 
investigated impact of late and early time cosmic models on the collapse of 
self-gravitating systems with metric as well as Palatini f{R) theory (Sharif 
and Yousaf 2013a,b,c,d, 2014a,c). 

A great deal of effort has been devoted to study the stability of stel¬ 
lar systems upon fluctuations. Galli and Koshelev (2011) studied a class of 
late-time cosmic evolution models with perturbations induced by inhomo¬ 
geneous energy density. Pinheiro and Ghan (2011) examined non-adiabatic 
anisotropic collapse accompanied by inhomogeneous density conhguration 
with and without shearing motion. Sharma and Tikekar (2012) explored 
shear-free spherical collapse with dissipation through heat to investigate the 
inhomogeneity effects during evolution. Sharif and his collaborators (Sharif 
and Yousaf 2012b,c, Sharif and Bhatti 2012a,b, 2014a,c, Sharif and Tahir 
2013) studied spherical, cylindrical and planar celestial models and analyzed 
the role of energy density inhomogeneity in the evolution of fluid parameters 
that characterize gravitational collapse. 

Herrera et al. (2009) explored spherical relativistic fluid conhgurations 
through scalar functions, i.e., Yt, Xt, Y^f and Xtf ■ Herrera et al. (2011) 
extended their results by invoking cosmological constant to examine the evo¬ 
lution of shear tensor and expansion scalar. They identihed Xtf as a factor 
describing inhomogeneity in the energy density. Herrera (2011) discussed 
stability of inhomogeneous density in anisotropic spherical fluid conhgura¬ 
tion with diffusion and free-streaming approximations. Recently, we have 
studied dynamics of spherical matter distribution with structure scalars and 
eR^ cosmology (Sharif and Yousaf 2014b). 

This paper investigates the role of polynomial f{R) gravity on the stabil¬ 
ity of homogeneous energy density with anisotropic and dissipative spherical 
matter. The paper is planned as follows. In section 2 , we discuss f{R) 
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formulation and relate matter variables with the Weyl scalar. Section 3 is 
devoted to construct scalar functions with a well-consistent polynomial f{R) 
model to obtain conservation and Ellis equations. In section 4, we consider 
various aspects of matter distribution to analyze density inhomogeneity. In 
the last section, we conclude our results. 


2 f{R) Formalism 


The gravitational part of the Einstein-Hilbert action in f{R) gravity is 

^ ^ j (1) 

where k and f{R) are coupling constant and a non-linear generic Ricci scalar 
function, respectively. The usual GR action can be retrieved by taking 
f{R) = R- The field equations, in metric formalism, are calculated by vary¬ 
ing Eq.(IT]) with respect to gap as follows 

Rapfu (V P 9aP^') fR '^Qapf iR^aPi (2) 

where and □ are the covariant derivative and d’Alembert operator, re¬ 
spectively. Equation ([2]), after some manipulations, can be expressed as 


K 

Gap = -^{TaP + Tap), (3) 

JR 

where 

Tap = I {VaVpfR + (/ - RIr)^ - DfRgap] , 

is the stress energy tensor which indicates gravitational contribution due to 

{D) 

f{R) terms. Under GR limit, i.e., f{R) —)■ R, Tap disappears identically. The 
system under consideration is modeled as a sphere with non-static spacetime 

ds^_ = r)dR — r)dr‘^ — C^{t, r){d6‘^ + sin^ 6d(j)^), (4) 


consisting of locally anisotropic pressure, dissipating in the diffusion (heat) 
and free streaming (null radiation) approximations. The corresponding stress- 
energy tensor is 


Tap — (/^ + P±)VaVp + {Pr ~ Pl.)XaXp ~ P±9oiP + da^P + Ra^P + VaQp, (5) 
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where /i, Pj_, P^, <?/3 and £, are the energy density, tangential and radial 
pressures, heat conducting vector and radiation density, respectively. More¬ 
over, and are the null four-vector, fluid four-velocity and radial 

unit four-vector, respectively. These quantities = ^(5q, x^ = = 

^^0 + ~ ill comoving coordinates obey 

1/“14 = -1, x“Xa = l, x“K = 0, 

= 0, rVa = -1, ria = 0. 


The expansion and shear scalars for Eq.([T]) are given by 

[2C b\ 

QA — I —-—h TT 1 , O'A — — 



( 6 ) 


where dot stands differentiation with respect to t. 
The metric /(P) field equations turn out to be 
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where prime represents differentiation with respect to r . The Misner-Sharp 
mass fnnction is given by (Misner and Sharp 1964) 


m 




The radial and proper derivative operators are defined respectively as follows 

( 12 ) 

The proper time rate of change of areal radius of the spherical system is 

(13) 

In terms of collapsing fluid velocity, Eq. fllip can be written as 

2m(t, r)"' 


Dc = —— Dt = —— 
^ C'dr ^ Adt' 


C 

U = DtC = — < 1 (for collapsing bodies). 


F-^' - 
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The time and radial mass variations can be followed from Eqs. 
flTTl) and flT^ as 
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where Pr = Pr + £, q = q + e and fi = jiP e. Integration of Eq. flT6D provides 


3m 

3k 

nr 




Jo 

fR 1 


{D) 

Too 

A2 


+ Q 


(D) 

Toi \ U 


AB E 




dr, (17) 


6 

















thereby relating mass function and other fluid variables with f{R) terms. 
The electric component of the Weyl tensor in terms of Xa and unit four 
velocity is given by 


where 


EaP — £ 


X.aX.j3 2 i,9ap T 


£ = 


which after using Eqs. 
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T 
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(mi) can be expressed as 


{D) (£>)' 

^ I - £. 


52 


C2 


(18) 


(19) 


where fl = Pr — P_l. This equation peculiarly relates mass function with 
Weyl scalar and all the fluid variables in f{P) gravity. 


3 Structure Scalars and Ellis Equations 

In this section, we hrst discuss a viable f{P) model and then construct 
structure scalars. We also write down conservation laws from the usual as 
well as effective stress energy tensors and develop the so called Ellis equations. 
We take a polynomial inflationary model given as follows (Huang 2014) 

/(/?) = + + (20) 

where e = and A„ is a dimension-free coupling parameter with n > 2. 
Here energy scale M is rehned in order to make unit normalization to the 
higher coefficient of P^ term. This model corresponds to the model with 
f{P) = P + P^/— 1 under A„ S> 1 while A„ —)■ 0 provides Starobinsky 
model (1980). In the limit A„ 1, i?” terms serve as a small correction to 
the inflationary P -|- eP^ model which of course makes the model expansion 
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around the Starobinsky model. It is interesting to mention here that the 
inflation induced by i? + gravity provides much different platform than 
that oi R + B? gravity and is close to topological inflation (Saidov 2010). All 
GR solutions can be found by taking limit f{R) —)■ R. 

To formulate f{R) structure scalars, we orthogonally split the Riemann 
tensor and propose tensors and Yap as (Herrera et ah 2011) 

where Rap^^i, = ap- These equations in terms of trace and trace-less 

components are given by 


^ap — -^Xxhap + Xtf 
Yap = -Yrhap + Ytf | 



( 21 ) 

( 22 ) 


We use Eqs. (ED, (ED, (HOD and fl2ni) - fl2^ with some manipulations to obtain 
the following scalar structures 
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(24) 
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(26) 

(26) 


where and (pp^ are given in Appendix A. It is well established in GR 

as well as in j{K) gravity that, one of the structure scalars Xt describes mat¬ 
ter energy density while its inhomogeneity is discussed with the help of Xtf 
only if the system evolves adiabatically alongwith e = 0. The scalar func¬ 
tions Ytf cuid Ip incorporating e terms control the evolutionary mechanisms 
of shearing and expansion rates of the system. 

The two independent components of the contracted Bianchi identities are 












which yield 

^ /A -N-S 2AqC' _ ^ . , , 

h + + (-Pr- + /i)— H——h 2(Pj_ + fx)— + Do{t,r) = 0, (28) 

Api A' (c b\ ( AC^' 

^ + ^+(A+/i)- + 2U + -jg + 2n^ + D,{t, r) = 0, (29) 

where Dq and Di are f{R) dark source terms given in Appendix A. Now 
we hnd two very important differential equations which play a pivotal in 
the stability analysis of inhomogeneous energy density. These two equations 
were firstly calculated by Ellis (2009) and then by Herrera et al. (2004) in 
GR. These equations are obtained by using Eqs.(j7 |) - ffT0|) . (IT^ . flTOD and (l20D 
as 
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(31) 


where (fq is mentioned in Appendix A. Both of the above equations reduce 
to GR (Herrera 2011) under e —)■ 0. 


4 Stability of Homogeneous Energy Density 

In this section, we discuss different factors affecting energy density homo¬ 
geneity in matter distribution with f{R) framework for different cases. We 
confine ourselves with present valued cosmological Ricci scalar, i.e., R = R. 
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4.1 Non-dissipative Fluids 

In this subsection, we perform our analysis with non-dissipative matter distri¬ 
bution with polynomial f{R) gravity model for dust, isotropic and anisotropic 
fluid configurations. 


4.1.1 Dust Cloud 


Here we take non-dissipative dust fluid with its geodesic motion which gives 
q = P± = Pr = 0 and A = 1. In this context, Eqs. fl30|) and give 
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By making use of Eqs.( 


and flB3P - flB6p in Eqs. (l3^ and 
—2KeAafiR 


c 4(R{1 + 2efi) + A„(2€fl)" ’ 




2KeRfi' 


C 4ei?(l + 2eR) + A„(2ei?)» 


(33) 
we obtain 

(34) 

(35) 


Equation fl3411 yields 


£ = 


-2KeR f* Aa^C^dt 


[4:eR{l + 2eR) + Xn{2eRY]C^' 


(36) 


which provides condition for the existence of homogeneity in the dust fluid. 
This states that non-dissipative homogeneous spherical matter configuration 
exists only if the system is conformally flat. Similarly, we can identify the 
Weyl scalar as an inhomogeneity factor from Eq. fl35p . Equation fl36l) also 
asserts that conformal flatness exists if the system evolves with vanishing 
shear scalar. 


10 





















4.1.2 Isotropic Fluid 

Now, we consider adiabatic spherical system with locally isotropic pressure. 
Under this scenario, Eqs. fl30|) and fl^ become 
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It is seen that Eq. fl38|l turns out to be same as Eq. fl33l) . thus showing that 
= 0 if and only if /i' = 0. Equation flTTI) after using Eqs.([HD fl25]l 
provides the following differential equation 


^ ^ —2KeAa{^ + P)R 

^ C ~ AeR{l + 2eR) + A„(2e.R)^ ’ 


(39) 


whose solution is 

^ ^ -2KeRjlAa{ii + P)C^dt 
[4e.R(l + 2eR) + \n{2eRY]C^' 

This argues that energy density of the system will be homogeneous as long 
as the system embodies shear-free motion. Thus the condition of locally 
isotropic pressure in the matter conhguration with constant Ricci scalar f{R) 
model does not affect the stability of homogeneous energy density found in 
the above case. Let us assume shear-free fluid so that Eq. (l3^ provides 


^ _ yy 

where a;(r) is an integration function. If the system is homogeneous initially, 
i.e., £{0,r) = 0, then ca = 0 yields £{t,r) = 0. Thus the above condition for 
homogeneous energy density will be valid from f = 0 to onward. However, 
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if the fluid is expanding such that £ has very small (non-zero) value at the 
initial stage, then it will remain as it is in all the evolutionary phases. If 
instead the system is contracting, then the Weyl scalar does not vanish for 
all time. 


4.1.3 Anisotropic Fluid 


This case corresponds to anisotropic but non-dissipating matter distribution, 
i.e. n 7 ^ 0 and g = 0. In this framework, Eqs. (l30|) and flMl) provide 
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We use Eqs.(l6]) and (128|1 in Eqs. (1TT]l and (1421) to obtain the following set of 
evolutionary equations 
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These equations in terms of structure scalar fl24l) can be written as 
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whose general solutions can be found respectively as 
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These equations indicate that quantity incorporating stability of inhomoge¬ 
neous energy density is one of the f{R) structure scalars, i.e., Xtf- Equation 
fl44|) shows that /x' = 0 if and only if Xtf vanishes, thereby showing Xtf as 
a factor of controlling inhomogeneity in anisotropic spherical system which is 
well-consistent with (Herrera et ah 2009, Herrera et ah 2011, Herrera 2011). 
Thus the inclusion of dark matter/energy effects in the evolving system do 
not disrupt the importance of Xtf- Also, the above expressions reduce to GR 
under the limit e —)■ 0. However, Eq. fH^ asserts that anisotropic pressure, 
f{R) model and shear scalar are responsible for the emergence of inhomoge¬ 
neous energy density in the matter distribution. 


4.2 Dissipative Dust Cloud 

To see the effects of radiation density and heat conducting vector in the 
inhomogeneous energy density, we assume geodesic case, i.e., = -Rl = 0 

with A = 1. Many authors (Herrera et ah 2004, Herrera 2011, Kolassis et ah 
1988, Govender et ah 1998, Thirukkanesh and Maharaj 2009, Naidu et ah 
2006) discussed spherical dissipative collapsing dust models with geodesics in 
order to explore dissipation effects through the system. In this case, Eqs. fl30l) 
and fin]) yield 
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It is found that /i' = 0 if and only if 4) = 0, indicating that <I> is responsible for 
fluid density inhomogeneity in the dust spherical system with free streaming 
and diffusion approximations. We use Eqs.dH]) and fl25)l in Eq. fl4^ to obtain 
<h evolution equation as follows 


_ 2_^ _ 

~ AeR{l + 2eR) + Xn{2eR)^ \BC BJ C 
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(49) 


This indicates that various fluid parameters affect the evolution of <I> in the 
self-gravitating system. We also see from the above relation that existence 
of inhomogeneous density depends upon two factors, i.e., dissipative param¬ 
eters and shear scalar. This describes that shearing scalar, radiation density 
and heat dissipation hold fundamental importance in the the study of inho¬ 
mogeneous matter distribution leading to gravitational collapse. 


5 Summary and Discussion 

This work analyzes various factors producing inhomogeneity in the energy 
density of the spherical self-gravitating celestial body in f{R) gravity. We 
have constructed structure scalars by orthogonally splitting the Riemann cur¬ 
vature tensor to obtain evolution equations using a viable inflationary f{R) 
model. We have discussed our analysis for non-dissipative dust, isotropic as 
well as anisotropic fluid conhgurations and dust cloud dissipating fluid. The 
results are concluded as follows. 
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• For non-dissipative dust and locally isotropic ideal matter cloud, it is 
seen from Eqs. fl35l) and fl38p that the system will encapsulate homoge¬ 
neous energy density if and only if the system is conformally flat. The 
extra f{R) degrees of freedom terms turn down contribution of S, thus 
relaxing conformal flatness condition. 

• In an adiabatic anisotropic spherical system, the density inhomogeneity 

is described in terms of pressure anisotropy which in turn controlled by 
one of the structure scalars, Xtf as mentioned in Eq. fl2T|) . Equation 
fl44p also establishes Xtf as an element of governing inhomogeneity 
in the system. This result is well-consistent with (Sharif and Yousaf 
2014b) under —)• 0 and (Herrera 2011) under e ^ 0 which correspond 

to solutions in i? -|- eR^ gravity and GR, respectively. 

• The quantity <I> is explored and identihed to be responsible for the 
emergence of inhomogeneity in energy density for geodesic radiating 
dust fluid. Extra curvature f{R) terms, dissipation parameters and 
shear scalar affect evolution of <F as described by Eq. fH^ . 

• All these results correspond to GR under e —>■ 0 (Herrera 2011). It is 
worth stressing that structure scalars obtained in Eqs. fl2^ - fl26P hold 
fundamental importance in the study of self-gravitating system. For 
A„ —>■ 0, scalar functions reduces for f{R) = R + eR^ cosmology (Sharif 
and Yousaf 2014b) while e —)■ 0 provides GR results (Herrera et ah 
2011 ). 
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The higher curvature terms Do and Di of Eqs. fl28P and (|29P are given as 
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